home *** CD-ROM | disk | FTP | other *** search
/ MacHack 2000 / MacHack 2000.toast / pc / The Hacks / MacHacksBug / Python 1.5.2c1 / Extensions / Numerical / Lib / LinearAlgebra.py < prev    next >
Encoding:
Text File  |  2000-06-23  |  11.0 KB  |  354 lines
  1. # This module is a lite version of LinAlg.py module which contains
  2. # high-level Python interface to the LAPACK library.  The lite version
  3. # only accesses the following LAPACK functions: dgesv, zgesv, dgeev,
  4. # zgeev, dgesvd, zdesvd, dgelss, zgelss.
  5.  
  6. import Numeric
  7. import copy
  8. import lapack_lite
  9.  
  10. # Error object
  11. LinAlgError = 'LinearAlgebraError'
  12.  
  13. # Helper routines
  14. _lapack_type = {'f': 0, 'd': 1, 'F': 2, 'D': 3}
  15. _lapack_letter = ['s', 'd', 'c', 'z']
  16. _array_kind = {'i':0, 'l': 0, 'f': 0, 'd': 0, 'F': 1, 'D': 1}
  17. _array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
  18. _array_type = [['f', 'd'], ['F', 'D']]
  19.  
  20. def _commonType(*arrays):
  21.     kind = 0
  22. #    precision = 0
  23. #   force higher precision in lite version
  24.     precision = 1
  25.     for a in arrays:
  26.         t = a.typecode()
  27.         kind = max(kind, _array_kind[t])
  28.         precision = max(precision, _array_precision[t])
  29.     return _array_type[kind][precision]
  30.  
  31. def _castCopyAndTranspose(type, *arrays):
  32.     cast_arrays = ()
  33.     for a in arrays:
  34.         if a.typecode() == type:
  35.             cast_arrays = cast_arrays + (copy.copy(Numeric.transpose(a)),)
  36.         else:
  37.             cast_arrays = cast_arrays + (copy.copy(
  38.                                        Numeric.transpose(a).astype(type)),)
  39.     if len(cast_arrays) == 1:
  40.             return cast_arrays[0]
  41.     else:
  42.         return cast_arrays
  43.  
  44. def _assertRank2(*arrays):
  45.     for a in arrays:
  46.         if len(a.shape) != 2:
  47.             raise LinAlgError, 'Array must be two-dimensional'
  48.  
  49. def _assertSquareness(*arrays):
  50.     for a in arrays:
  51.         if max(a.shape) != min(a.shape):
  52.             raise LinAlgError, 'Array must be square'
  53.     
  54.  
  55. # Linear equations
  56.         
  57. def solve_linear_equations(a, b):
  58.     one_eq = len(b.shape) == 1
  59.     if one_eq:
  60.         b = b[:, Numeric.NewAxis]
  61.     _assertRank2(a, b)
  62.     _assertSquareness(a)
  63.     n_eq = a.shape[0]
  64.     n_rhs = b.shape[1]
  65.     if n_eq != b.shape[0]:
  66.         raise LinAlgError, 'Incompatible dimensions'
  67.     t =_commonType(a, b)
  68. #    lapack_routine = _findLapackRoutine('gesv', t)
  69.     if _array_kind[t] == 1: # Complex routines take different arguments
  70.         lapack_routine = lapack_lite.zgesv
  71.     else:
  72.         lapack_routine = lapack_lite.dgesv
  73.     a, b = _castCopyAndTranspose(t, a, b)
  74.     pivots = Numeric.zeros(n_eq, 'l')
  75.     results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
  76.     if results['info'] > 0:
  77.         raise LinAlgError, 'Singular matrix'
  78.     if one_eq:
  79.         return copy.copy(Numeric.ravel(b))
  80.     else:
  81.         return copy.copy(Numeric.transpose(b))
  82.  
  83.  
  84. # Matrix inversion
  85.  
  86. def inverse(a):
  87.     return solve_linear_equations(a, Numeric.identity(a.shape[0]))
  88.  
  89. # Eigenvalues
  90.  
  91. def eigenvalues(a):
  92.     _assertRank2(a)
  93.     _assertSquareness(a)
  94.     t =_commonType(a)
  95.     real_t = _array_type[0][_array_precision[t]]
  96.     a = _castCopyAndTranspose(t, a)
  97.     n = a.shape[0]
  98.     dummy = Numeric.zeros((1,), t)
  99.     lwork = max(1,6*n) # minimum value: max(1,3*n) real, max(1,2*n) complex
  100.     work = Numeric.zeros((lwork,), t)
  101.     if _array_kind[t] == 1: # Complex routines take different arguments
  102.         lapack_routine = lapack_lite.zgeev
  103.         w = Numeric.zeros((n,), t)
  104.         rwork = Numeric.zeros((n,),real_t)
  105.         results = lapack_routine('N', 'N', n, a, n, w,
  106.                                  dummy, 1, dummy, 1, work, lwork, rwork, 0)
  107.     else:
  108.         lapack_routine = lapack_lite.dgeev
  109.         wr = Numeric.zeros((n,), t)
  110.         wi = Numeric.zeros((n,), t)
  111.         results = lapack_routine('N', 'N', n, a, n, wr, wi,
  112.                                  dummy, 1, dummy, 1, work, lwork, 0)
  113.         if Numeric.logical_and.reduce(Numeric.equal(wi, 0.)):
  114.             w = wr
  115.         else:
  116.             w = wr+1j*wi
  117.     if results['info'] > 0:
  118.         raise LinAlgError, 'Eigenvalues did not converge'
  119.     return w
  120.  
  121.  
  122. # Eigenvectors
  123.  
  124. def eigenvectors(a):
  125.     """eigenvectors(a) returns u,v  where u is the eigenvalues and
  126. v is a matrix of eigenvectors with vector v[i] corresponds to 
  127. eigenvalue u[i].  Satisfies the equation matrixmultiply(a, v[i]) = u[i]*v[i]
  128. """
  129.     _assertRank2(a)
  130.     _assertSquareness(a)
  131.     t =_commonType(a)
  132.     real_t = _array_type[0][_array_precision[t]]
  133.     a = _castCopyAndTranspose(t, a)
  134.     n = a.shape[0]
  135.     lwork = max(1,8*n) # minimum value: max(1,4*n) real, max(1,2*n) complex
  136.     work = Numeric.zeros((lwork,), t)
  137.     dummy = Numeric.zeros((1,), t)
  138.     if _array_kind[t] == 1: # Complex routines take different arguments
  139.         lapack_routine = lapack_lite.zgeev
  140.         w = Numeric.zeros((n,), t)
  141.         rwork = Numeric.zeros((lwork,),real_t)
  142.         v = Numeric.zeros((n,n), t)
  143.         results = lapack_routine('N', 'V', n, a, n, w,
  144.                                   dummy, 1, v, n, work, lwork, rwork, 0)
  145.     else:
  146.         lapack_routine = lapack_lite.dgeev
  147.         wr = Numeric.zeros((n,), t)
  148.         wi = Numeric.zeros((n,), t)
  149.         vr = Numeric.zeros((n,n), t)
  150.         results = lapack_routine('N', 'V', n, a, n, wr, wi,
  151.                                   dummy, 1, vr, n, work, lwork, 0)
  152.         if Numeric.logical_and.reduce(Numeric.equal(wi, 0.)):
  153.             w = wr
  154.             v = vr
  155.         else:
  156.             w = wr+1j*wi
  157.             v = Numeric.array(vr,Numeric.Complex)
  158.             ind = Numeric.nonzero(
  159.                           Numeric.equal(
  160.                               Numeric.equal(wi,0.0) # true for real e-vals
  161.                                        ,0)          # true for complex e-vals
  162.                                  )                  # indices of complex e-vals
  163.             for i in range(len(ind)/2):
  164.                 v[ind[2*i]] = vr[ind[2*i]] + 1j*vr[ind[2*i+1]]
  165.                 v[ind[2*i+1]] = vr[ind[2*i]] - 1j*vr[ind[2*i+1]]
  166.     if results['info'] > 0:
  167.         raise LinAlgError, 'Eigenvalues did not converge'
  168.     return w,v
  169.  
  170.  
  171.    
  172. # Singular value decomposition
  173.  
  174. def singular_value_decomposition(a, full_matrices = 0):
  175.     _assertRank2(a)
  176.     n = a.shape[1]
  177.     m = a.shape[0]
  178.     t =_commonType(a)
  179.     real_t = _array_type[0][_array_precision[t]]
  180.     a = _castCopyAndTranspose(t, a)
  181.     if full_matrices:
  182.     nu = m
  183.     nvt = n
  184.     option = 'A'
  185.     else:
  186.     nu = min(n,m)
  187.     nvt = min(n,m)
  188.     option = 'S'
  189.     s = Numeric.zeros((min(n,m),), real_t)
  190.     u = Numeric.zeros((nu, m), t)
  191.     lwork = 10*max(m,n) # minimum value: max(3*min(m,n)+max(m,n),5*min(m,n)-4)
  192.     work = Numeric.zeros((lwork,), t)
  193.     vt = Numeric.zeros((n, nvt), t)
  194.     if _array_kind[t] == 1: # Complex routines take different arguments
  195.         lapack_routine = lapack_lite.zgesvd
  196.         rwork = Numeric.zeros((max(3*min(m,n),5*min(m,n)-4),), real_t)
  197.         results = lapack_routine(option, option, m, n, a, m, s, u, m, vt, nvt,
  198.                                  work, lwork, rwork, 0)
  199.     else:
  200.         lapack_routine = lapack_lite.dgesvd
  201.         results = lapack_routine(option, option, m, n, a, m, s, u, m, vt, nvt,
  202.                                  work, lwork, 0)
  203.     if results['info'] > 0:
  204.         raise LinAlgError, 'SVD did not converge'
  205.     return copy.copy(Numeric.transpose(u)), s, copy.copy(Numeric.transpose(vt))
  206.  
  207.  
  208. # Generalized inverse
  209.  
  210. def generalized_inverse(a, rcond = 1.e-10):
  211.     u, s, vt = singular_value_decomposition(a, 0)
  212.     m = u.shape[0]
  213.     n = vt.shape[1]
  214.     cutoff = rcond*Numeric.maximum.reduce(s)
  215.     for i in range(min(n,m)):
  216.         if s[i] > cutoff:
  217.             s[i] = 1./s[i]
  218.     else:
  219.         s[i] = 0.;
  220.     return Numeric.dot(Numeric.transpose(vt),
  221.                s[:, Numeric.NewAxis]*Numeric.transpose(u))
  222.  
  223. # Determinant
  224.  
  225. def determinant(a):
  226.     _assertRank2(a)
  227.     _assertSquareness(a)
  228.     t =_commonType(a)
  229.     a = _castCopyAndTranspose(t, a)
  230.     n = a.shape[0]
  231.     if _array_kind[t] == 1:
  232.     lapack_routine = lapack_lite.zgetrf
  233.     else:
  234.     lapack_routine = lapack_lite.dgetrf
  235.     pivots = Numeric.zeros((n,), 'l')
  236.     results = lapack_routine(n, n, a, n, pivots, 0)
  237.     sign = Numeric.add.reduce(Numeric.not_equal(pivots,
  238.                         Numeric.arrayrange(1, n+1))) % 2
  239.     return (1.-2.*sign)*Numeric.multiply.reduce(Numeric.diagonal(a))
  240.  
  241.  
  242. # Linear Least Squares 
  243.         
  244. def linear_least_squares(a, b, rcond=1.e-10):
  245.     """solveLinearLeastSquares(a,b) returns x,resids,rank,s 
  246. where x minimizes 2-norm(|b - Ax|) 
  247.       resids is the sum square residuals
  248.       rank is the rank of A
  249.       s is an rank of the singual values of A in desending order
  250.  
  251. If b is a matrix then x is also a matrix with corresponding columns.
  252. If the rank of A is less than the number of columns of A or greater than
  253. the numer of rows, then residuals will be returned as an empty array
  254. otherwise resids = sum((b-dot(A,x)**2).
  255. Singular values less than s[0]*rcond are treated as zero.
  256. """
  257.     one_eq = len(b.shape) == 1
  258.     if one_eq:
  259.         b = b[:, Numeric.NewAxis]
  260.     _assertRank2(a, b)
  261.     m  = a.shape[0]
  262.     n  = a.shape[1]
  263.     n_rhs = b.shape[1]
  264.     ldb = max(n,m)
  265.     if m != b.shape[0]:
  266.         raise LinAlgError, 'Incompatible dimensions'
  267.     t =_commonType(a, b)
  268.     real_t = _array_type[0][_array_precision[t]]
  269.     bstar = Numeric.zeros((ldb,n_rhs),t)
  270.     bstar[:b.shape[0],:n_rhs] = copy.copy(b)
  271.     a,bstar = _castCopyAndTranspose(t, a, bstar)
  272.     lwork = 8*min(n,m) + max([2*min(m,n),max(m,n),n_rhs]) 
  273.     s = Numeric.zeros((min(m,n),),real_t)
  274.     work = Numeric.zeros((lwork,), t)
  275.     if _array_kind[t] == 1: # Complex routines take different arguments
  276.         lapack_routine = lapack_lite.zgelss
  277.         rwork = Numeric.zeros((5*min(m,n)-1,), real_t)
  278.         results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
  279.                         0,work,lwork,rwork,0 )
  280.     else:
  281.         lapack_routine = lapack_lite.dgelss
  282.         results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
  283.                         0,work,lwork,0 )
  284.     if results['info'] > 0:
  285.         raise LinAlgError, 'SVD did not converge in Linear Least Squares'
  286.     resids = Numeric.array([],t)
  287.     if one_eq:
  288.         x = copy.copy(Numeric.ravel(bstar)[:n])
  289.         if (results['rank']==n) and (m>n):
  290.             resids = Numeric.array([Numeric.sum((Numeric.ravel(bstar)[n:])**2)])
  291.     else:
  292.         x = copy.copy(Numeric.transpose(bstar)[:n,:])
  293.         if (results['rank']==n) and (m>n):
  294.             resids = copy.copy(Numeric.sum((Numeric.transpose(bstar)[n:,:])**2))
  295.     return x,resids,results['rank'],copy.copy(s[:min(n,m)]) 
  296.  
  297.  
  298. if __name__ == '__main__':
  299.     from Numeric import *
  300.  
  301.     def test(a, b):
  302.  
  303.         print "All numbers printed should be (almost) zero:"
  304.  
  305.         x = solve_linear_equations(a, b)
  306.         check = b - matrixmultiply(a, x)
  307.         print check
  308.  
  309.  
  310.         a_inv = inverse(a)
  311.         check = matrixmultiply(a, a_inv)-identity(a.shape[0])
  312.         print check
  313.  
  314.  
  315.         ev = eigenvalues(a)
  316.  
  317.         evalues, evectors = eigenvectors(a)
  318.         check = ev-evalues
  319.         print check
  320.  
  321.         evectors = transpose(evectors)
  322.         check = matrixmultiply(a, evectors)-evectors*evalues
  323.         print check
  324.  
  325.  
  326.         u, s, vt = singular_value_decomposition(a)
  327.         check = a - Numeric.matrixmultiply(u*s, vt)
  328.         print check
  329.  
  330.  
  331.         a_ginv = generalized_inverse(a)
  332.         check = matrixmultiply(a, a_ginv)-identity(a.shape[0])
  333.         print check
  334.  
  335.  
  336.         det = determinant(a)
  337.         check = det-multiply.reduce(evalues)
  338.         print check
  339.  
  340.         x, residuals, rank, sv = linear_least_squares(a, b)
  341.         check = b - matrixmultiply(a, x)
  342.         print check
  343.         print rank-a.shape[0]
  344.         print sv-s
  345.  
  346.     a = array([[1.,2.], [3.,4.]])
  347.     b = array([2., 1.])
  348.     test(a, b)
  349.  
  350.     a = a+0j
  351.     b = b+0j
  352.     test(a, b)
  353.  
  354.